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ABSTRACT 

Circumbinary planetary systems recently discovered by Kepler represent an important testbed for 
planet formation theories. Planetesimal growth in disks around binaries has been expected to be 
inhibited interior to ~ 10 AU by secular excitation of high relative velocities between planetesimals, 
leading to their collisional destruction (rather than agglomeration). Here we show that gravity of 
the gaseous circumbinary disk in which planets form drives fast precession of both the planetesimal 
and binary orbits, resulting in strong suppression of planetesimal eccentricities beyond 2-3 AU and 
making possible growth of 1 — 10 2 km objects in this region. The precise location of the boundary 
of accretion-friendly region depends on the size of the inner disk cavity cleared by the binary torques 
and on the disk mass (even 0.01 M disk strongly suppresses planetesimal excitation), among other 
things. Precession of the orbit of the central binary, enhanced by the mass concentration naturally 
present at the inner edge of a circumbinary disk, plays key role in this suppression, which is a feature 
specific to the circumbinary planet formation. 

Subject headings: planets and satellites: formation — protoplanetary disks — planetary systems — 
binaries: close 



1. INTRODUCTION. 

One of the most intriguing findings of the Kepler mis- 
sion is the discovery of circumbinary planets around a 
sample of close binaries, affectionately termed Tatooines 
(Doyle et al. 2011; Welsh et al. 2012; Orosz et al. 2012a,b; 
Schwamb et al. 2012). Six such systems are known at 
the moment, with one of them, Kepler-47, harboring two 
planets (see Table [1]). At least some of the Kepler cir- 
cumbinary planets are likely in the Saturn or Jupiter 
mass regime. Semi-major axes of their orbits are typi- 
cally small, a p i < 1 AU, and close to the limit of dynam- 
ical stability (Holman & Wiegert 1999). 

These systems provide interesting targets for testing 
our understanding of planet formation. It is generally 
agreed that in-situ formation of such planets is impos- 
sible because of the strong dynamical excitation due 
to central binary at their present locations (Meschiari 
2012a; Paardekooper et al. 2012). The question typically 
addressed is at what separation could these planets form, 
subsequent to which they have migrated in. The bottle- 
neck for planet formation here is the growth of 1 — 10 2 
km planetesimals, which is impossible if planetesimals 
have high eccentricities — instead of merging they get de- 
stroyed in high-speed collisions. Since the excitation by 
the binary is a rather long-range effect (in conventional 
secular theory eccentricity decays with distance only as 
see equation ([23])). it might easily be the case that 
the conditions for planetesimal growth are realized only 
beyond 10 AU (Moriwaki & Nakagawa 2004, hereafter 
MN04; Scholl et al. 2007). However, at such separa- 
tions the timescale for growing massive cores capable of 
forming giant planets by core accretion may become pro- 
hibitive, especially if the lifetimes of circumbinary disks 
are shorter than around single stars. The latter seems to 
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be the case at least for the circumstellar disks in binaries 
(Cieza et al. 2009; Kraus et al. 2012). 

Different ideas were proposed to alleviate this plan- 
etesimal fragmentation issue. Orbital alignment due to 
gas drag in presence of secular forcing has been sug- 
gested to facilitate growth of 1 — 10 km bodies (Scholl et 
al. 2007), but this effect is size-dependent and does not 
work well for a broad spectrum of planetesimal masses. 
Paardekooper et al. (2012) explored very efficient accre- 
tion of dust by growing planetesimals as another way of 
bypassing the fragmentation problem. 

In this work we show that gravitational effect of the 
circumbinary disk, in which planetesimals are immersed, 
modifies their secular excitation and significantly lowers 
their eccentricities beyond 2-4 AU from the star. Similar 
mechanism was invoked in Rafikov (2012b) to explain the 
origin of circumstellar planets on wide orbits (~ 2 AU) in 
small separation binaries, ab ~ 20 AU. A unique feature 
of the circumbinary planet formation is that the effect 
of the disk gravity on the central binary often turns out 
being more important than the direct effect of the disk on 
planetesimal orbits. We now describe this idea in more 
detail. 

2. SECULAR EVOLUTION. 

We consider motion of massless planetesimals around a 
central binary. Binary has components with masses M p 
and M s < M p (mass ratio \i = M s jM\, < 1, where Mf, = 
M p + M s is the binary mass), its semimajor axis and 
eccentricity are ab and eb correspondingly. For simplicity 
we will assume the disk to be axisymmetric with respect 
to the binary barycenter, i.e. the disk surface density is 
E(r), where r is the distance from the barycenter. Since 
we are primarily interested in the effect of the binary on 
planetesimal dynamics we neglect gas drag in this work 
(see Planetesimals start on circular orbits and we 
are interested in their eccentricity evolution driven by 
the time-dependent and non-axisymmetric potential of 
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Table 1 

Circumbinary planetary systems 



System 


a b 
(AU) 




M p 
(Mq) 


M s 
(M ) 


a P i 
(AU) 


^form 

(AU) 


Kepler-16 1 


0.22 


0.16 


0.69 


0.2 


0.7 


2.3-4.4 


Kepler-34 2 


0.23 


0.52 


1.05 


1.02 


1.1 


2.1-4.3 


Kepler-35 2 


0.18 


0.14 


0.89 


0.81 


0.6 


1.7-3.5 


Kepler-38 3 


0.15 


0.1 


0.95 


0.25 


0.46 


1.9-3.5 


Kepler-47 4 


0.084 


0.02 


1.04 


0.36 


0.29 b 


1.0-1.9 


KIC 4862625 5 


0.17 


0.21 


1.38 


0.39 


0.63 


2.7-5.1 



1 Doyle et al. (2011); 2 Welsh et al. (2012); 3 Orosz et al. (2012b); 
4 Orosz et al. (2012a); 5 Schwamb et al. (2012) 
a Inner edge of accretion-friendly zone, see JS] 
b Semi-major axis of the inner planet 



the binary. 

Following MN04 we write down the equation of plan- 
etcsimal motion as 



d 2 r _ G{M p + M s 



dt 2 r 3 
where the disturbing function R is 
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Here Ud is the disk potential and r p and r s are the vectors 
to primary and secondary from the barycenter of the 
binary. 

Using Murray & Dermott (1999) we expand the dis- 
turbing function to second order in planetesimal and bi- 
nary eccentricities e and e b , retaining terms up to 0(e 2 ) 
and 0(ee b ), and then average it over the binary and plan- 
etesimal mean longitudes, thus eliminating short-period 
perturbations. We additionally expanded Laplace coef- 
ficients assuming a b /a <C 1, where a is the planetesimal 
semi-major axis. The resulting secular disturbing func- 
tion R sec is 
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where we dropped insignificant e-independent terms. 
Here n b = (GM b /a 3 b )^ 2 is the binary mean motion, 
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is the precession frequency of planetesimal orbit due to 
the disk potential, and n = {GM b / a 3 ) 1 / 2 is the mean 
motion around a point mass M p + M s . To lowest order 
in e these expressions coincide with secular expansion of 
MN04, if we set w d = 0. 



2.1. Evolution equations. 

We now follow standard procedure (Murray & Der- 
mott 1999) and introduce eccentricity vector e = (k, h) = 



(e cos to, e sin w). Defining 

3 ~ 2 
4 
15 
16 

we can write 
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(ft, 2 + k 2 ) + B(kcosw b + hsmw b ).(7) 



This expression is accurate to 0(e 2 ). Evolution equa- 
tions for h and k (Murray & Dermott 1999) attain a 
relatively simple form 
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In the disk-free case (wd = 0) one recovers the secular 
evolution equations from MN04. 

Now we introduce an important modification to the 
setup used in MN04. We assume that w b is not con- 
stant but linearly increases with time at constant rate 
Wd, which we specify later in £j3j i.e. w b — w^t. 
This makes forcing term in equations jSJ) and (j9|) time- 
dependent, but still permits analytical solution in the 
form e(t) = e froc (t) + e {olced (t), where 
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, ^ ■ (12) 

A + Wd — W b 

Thus, free eccentricity vector ef ree rotates at a rate A+Wd 
around the endpoint of vector ef orC cd , which itself rotates 
about the origin with the rate w b . Setting Wd = to& = 
brings us back to the MN04 solution. 

Planetesimals starting on circular orbits have ef rcc = 
eforced and reach the maximum eccentricity of 

2B . . 

e = 2ef orce d = , ; . — (13) 

A + w d — Wb 

in the course of their secular evolution. 

3. DISK MODEL. 

Structure of circumbinary disks is different from that 
of protoplanetary disks around single stars. As a result 
of viscous evolution the latter are expected to have mass 
accretion rate M independent of radius. In the case of 
circumbinary disk the torque due to binary stops the 
inward flow of matter and truncates the disk at inner 
radius Ti n , interior to which S and M are small. Simu- 
lations find that rj„ 2a b for binaries with mass ratio 
Ms/Mp ~ 1 (MacFadyen & Milosavljevic 2008). 

Injection of angular momentum at its inner edge causes 
the disk to evolve into a configuration, in which the vis- 
cous angular momentum flux 



Fj = 37ri/£fir 2 , 



(14) 
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rather than M, is constant with radius (Pringle 1991; 
Ivanov et al. 1999; Rafikov 2012a). Detailed descrip- 
tion of circumbinary protoplanetary disk properties will 
be provided elsewhere (Garmilla & Rafikov, in prepa- 
ration); for the purposes of this paper we will assume 
that the disk is passive, i.e. heated predominantly by 
the combined light of the binary components. Assuming 
a-model for the viscosity v one finds 



a, AU 



E = 



F, 



3nac 2 r 2 



(15) 



where c s is a sound speed. If a is independent of radius, 
and midplane disk temperature scales as T(r) oc r~ k , 
then the constant Fj disk has density profile 



E(r) 



p = 2-k, 
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where Ej„ = E(ri n ). Passive disks typically have k w 
1/2, in particular Chiang & Goldreich (1997) find k = 
3/7. For this reason we will take p = 3/2 in this work, 
which is similar to the E slope of the Minimum Mass 
Solar Nebula (Hayashi 1981), and is different from p w 1 
expected for a passive constant M disk (Rafikov 2012b). 

Most of the mass in a constant Fj disk is contained 
in its outer regions. Assuming that p = 3/2 profile is 
maintained all the way to the outer radius r Q we find 
that 
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where M d is the disk mass, r a 30 = r o /(30AU), r% — 
r/(lAU). 

We show in Appendix [A] that disk with density profile 
([I?]) gives rise to apsidal precession of planetesimal orbits 
at the rate 
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(where K d w 4.4), as well as apsidal precession of the 
central binary at the rate 
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where cj) ~ 0.5, see equation (|A3|) . 

Equation (|20p assumes disk to be sharply truncated at 
r.i n , while in reality E smoothly (but rapidly) goes to zero 
near Ti n (MacFadyen & Milosavljevic 2008), which lowers 
the amount of mass near the binary and may reduce Wb- 
To account for this we will sometimes consider wider disk 
cavity, e.g. ri n = 3a?, (Pelupessy & Portegies Zwart 2012) 
instead of the more conventional 2af,. 

4. PLANETESIMAL VELOCITIES. 

Equation (TT2l shows that planetesimal eccentricity is 
determined by A, w d and Wb- We plot the behav- 
ior of these frequencies as a function of r in Figure 
[Q Using equations ©, (JTSJ) and ([20]) we find that 
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Figure 1. Characteristic precession frequencies of the problem 
as a function of distance from the binary (in units of aj on the 
lower axis and in AU on the upper), for different values of the 
relative disk mass M^/Mi, and the size of the inner disk cavity 
rin/a-by as labeled on panels. Calculation assumes Mj, = M@, 
/i = 0.21, a b = 0.2 AU, and e b = 0.2. Shown are A/n (solid), 
\zu,l\/ n {dotted), and m^ln (dashed), marked in panel (a). Note 
that for small M^/Mf, and Tin/a^ region where \"A d \ dominates 
over other frequencies disappears. 



Wb/A oc (r/ab) 7 / 2 and Wbj\w d \ oc r/ab. Both ratios 
increase with r so that Wb ^ A for r > ta, where 
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with ab,o.2 = at/(0.2AU). Also, Wb > (rorfj for r > r^, 
where (setting ifd = 4.4 and = 0.5) 

5/2 

(22) 



Thus, for a particular set of parameters adopted in 
these estimates one finds that Wb dominates the behavior 
of planetesimal eccentricity beyond (10 — 20)a&. Depend- 
ing on M d /Mb and r,- n /a6 there may exist intermediate 
region around r ~ 10ah, where relative precession of the 
planetesimal and binary orbits is dominated by w d . Such 
region disappears for lower M d and smaller inner cavity 
size ri n , see Figure [TJo. 

In Figure [2] we illustrate radial dependence of the char- 
acteristic planetesimal eccentricity e(r) computed with 
equation (fTB"]) . One feature that is obvious in these plots 
is the secular resonance at ~ 1 — 2 AU, where A is equal 
to zub + \w d \ and our solution presented in §2.11 breaks 
down. This region of enhanced excitation is quite nar- 
row, and just outside of it e(r) rapidly declines with r. 

We also show in this Figure asymptotic behavior of 
eccentricity e A (r) found when the relative precession of 
planetesimal and binary orbits is dominated by the po- 
tential of the secondary (A > |d7rf|, vJb), see ([5]), ([6]): 
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Fig ure 2. Maximum planetesimal eccentricity e(r) (solid curve, 
Eq. 1131 ) as a function of radius, for different values of M^/Mf, and 
Tin/o-b- Dot-dashed and dashed lines illustrate asymptotic behav- 
ior of eccentricity given by e A (Eq. [23]) and e b (Eq. [24) ^ . All 
these curves scale linearly with the binary eccentricity e;,, assumed 
equal to 0.2 in this calculation (also Mf, = Mq, ji = 0.21, a b = 0.2 
AU). Dotted line is the eccentricity e co u (Eq. [28 j ) below which 
10 km objects (density p = 3 g cm -3 ) can grow in collisions, ac- 
cording to the criterion {27} ■ Planetesimal growth is unimpeded 
by fragmentation when e(r) < e co u(r) (region to the right from 
vertical bars with arrows). 
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When relative precession is dominated by the binary 
precession (w b > A, Itu,^), equation (|20|) predicts e — > 
e b = 2B/w b , where 
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Here = jtx(l — /i)(l — 2/z) and for < \i < 0.5 the 
maximum value of tp w 0.096 is achieved at /i w 0.21. 

Between these two limits, at r ~ 1 AU, disk-driven 
planetesimal precession may dominate (see Figure [1]), 
but usually marginally. For this reason we do not show 
asymptotic scaling for this regime. 

At large separations, of order several AU, e(r) — > e b (r) 
and falls off very steeply with r. At these separations the 
simple formula (|23[) not accounting for the gravitational 
effects of the disk overestimates planetesimal eccentricity 
by more than an order of magnitude. 

5. SHORT-PERIOD ECCENTRICITY VARIATIONS. 



Planetesimal motion is affected not only by the ex- 
plicitly time-independent, secular part of the disturbing 
function (|3]), but also by the short-term perturbations 
varying on timescales ~ and The former aver- 
age out to zero over the planetesimal orbital motion, but 
the latter were suggested to affect planetesimal dynam- 
ics. 

In particular, MN04 and Paardekooper et al. (2012) 
found that even for circular binaries (e b = 0), when 
secular excitation is absent (see equations [6] and [13]), 
time-dependent contributions to the disturbing function 
varying on the planetesimal orbital timescale n~ x (i.e. 
averaged over the binary period) still result in eccentric- 
ity evolution. This result is surprising since in e b — case 
the potential of the binary averaged over the fast binary 
orbital timescale is time-independent and axisymmetric. 
Consequently, both energy and angular momentum of 
the planetesimal must be conserved precluding its eccen- 
tricity evolution. 

The unexpected finding of MN04 and Paardekooper et 
al. (2012) is most likely related to their choice of a refer- 
ence circular orbit and osculating orbital elements. The 
former was defined as the circular orbit in a Keplerian 
potential for the mass M b = M p + M s . However, one can 
show that the true axisymmetric part of the binary (plus 
disk) potential to lowest order in a b /r is 



Um=o(r) = U d (r)- 



GM b 



1 



1 + 7^ 
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Circular orbits in this potential have higher circular 
speed than orbits in a Keplerian potential for mass M b at 
the same distance. Thus, a particle initialized on "circu- 
lar" orbit, assuming —GM b /r potential, is in fact started 
at the apoapse of eccentric orbit in the true potential 
(|26p . Not surprisingly, its orbit will remain eccentric, 
with eccentricity determined ultimately by the difference 
between the circular velocities in U m =o and the M b point 
mass potentials. It is trivial to show that this fictitious 
eccentricity is (3/4)/x(l — fi)(a b /r) 2 , which is in perfect 
agreement with the numerical calculations of MN04 and 
analytical result Paardekooper et al. (2012) for e b = 0. 
The importance of proper definition of osculating orbital 
elements has been previously emphasized by Marzari et 
al. (2008). 

Careful analysis in the e b ^ case shows that the dis- 
turbing function contains terms varying on planetesimal 
orbital timescale, which formally result in eccentricity of 
order — fj){a b /r) 2 e b . Even taken at face value, this 
eccentricity is most likely too small to affect the results 
in f}4] whether it produces significant relative velocity be- 
tween colliding objects is even less obvious. We leave the 
detailed study of the short-term eccentricity variations 
for the future. 

6. IMPLICATIONS FOR PLANET FORMATION. 

We now address the issue of fragmentation barrier for 
planetesimal growth. Based on work of Leinhardt & 
Stewart (2012) we estimate that a gravity-dominated 
body (e.g. a rubble pile) of radius d survives in a col- 
lision with a body of equal size whenever the collision 
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velocity v co u (at large separation) satisfies a simple con- 
dition (Rafikov 2012b) 



Vcoii < 2v e 



(27) 



where the escape speed from the surface of an object of 
radius d and bulk density p is v esc — [(8n/3)Gp] l / 2 d. 
This condition is similar to the one used in MN04. It 
can be translated into the constraint on the maximum 
planetesimal eccentricity e co u ~ 2v esc /vK (where Vr is 
the Keplerian speed), at which an equal- mass collision 
does not result in the net loss of mass: 
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where M fe>1 = M b /M Q , p 3 = p/(3 g cm" 3 ), and d w = 
d/(10km). 

We will now assume that if (1) the characteristic plan- 
etesimal eccentricity given by equation (|13j) is below 
&coii{T^d s ) at some distance r and (2) planetesimals are 
strength-dominated below the radius d s , then fragmen- 
tation barrier at this separation r is bypassed. 

To be specific, we take d s = 10 km in our study. 
This may seem somewhat large since collisionally assem- 
bled objects may be rubble piles. On the other hand, 
more sophisticated fragmentation criteria accounting for 
the size spectrum of colliding objects (i.e. incorporat- 
ing more than just equal-mass collisions) typically find 
our survival criterion (l27l) too restrictive (Thebault 2011; 
Rafikov 2012b). This may justify relatively large value 
of d s = 10 km. 

Figure [2] shows that with this value of d s fragmentation 
barrier does not get in the way of planetesimal growth 
at separations > 2 — 4 AU. Inner radius of the accretion- 
friendly zone af orm depends mainly on the size of the 
inner cavity in the disk and to some extent on the disk 
mass. A rough estimate of af or m (typically an overesti- 
mate by up to a factor of 2) can be obtained by equating 
e b and e coU : 
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Larger ri n /a b means less disk mass near the binary and 
slower binary precession. Larger Md accelerates disk- 
induced precession of both the planetesimal orbit and 
the binary. Larger value of critical planetesimal size d s 
also shifts af or m closer to the binary. 

In the absence of disk-driven precession, using e A in- 
stead of e b in the growth condition (|2"7|) to calculate Of orm , 
we would obtain af orm w 90a;, « 17 AU (crossing of dot- 
ted and dot-dashed lines in Figure [2]), pushing accretion- 
friendly zone much further from the binary, in agreement 
with MN04. Formation of cores massive enough to trig- 
ger core accretion (~ 10M®) prior to disk dispersal is 
more problematic at this separation than at 2-3 AU, be- 
cause of the longer dynamical timescale. Thus, by ex- 
tending inward the region, where planetesimals can grow 



effectively, the gravitational effect of the disk on plan- 
etesimal secular evolution facilitates circumbinary planet 
formation via the reduction of the planetary accretion 
timescale. 

In Table Q] we show af orm computed using parame- 
ters of actual Kepler circumbinary planetary systems. 
Minimum and maximum values shown correspond to 
M d /M b = 0.1, r m /a b = 2 and M d /M b = 0.01, r m /a b = 
3, correspondingly. In Kepler-47 systenfl secular excita- 
tion is additionally suppressed because e b = 0.02 <C 1 
(Orosz et al. 2012a). In Kepler-34 and Kepler-35 it is 
reduced because p is very close to 0.5: p — 0.493 and 
0.477, correspondingly (Welsh et al. 2012). Compared to 
the nominal p = 0.2 case, we find e A in these systems to 
be suppressed by « 40 and w 13, respectively, while e d 
and e b are lowered by ~ 27 and 9, correspondingly. This 
shifts accretion-friendly zone closer to the binary. 

Nevertheless, even accounting for the gravitational ef- 
fect of the disk, in-situ formation still does not seem 
viable for the Kepler circumbinary planets, which have 
Oform > CLpi (Table [T]). Thus, some form of inward mi- 
gration is still needed to deliver these planets to their 
current locations. 

We also note that gas drag is not important for re- 
solving the fragmentation barrier issue in Kepler cir- 
cumbinary systems: fast relative precession of planetes- 
imal and binary orbits (mainly due to w b S> A) makes 
planetesimal apsidal alignment inefficient, contrary to 
standard expectation without disk gravity (Scholl et al. 
2007). Calculations similar to the one in Rafikov (2012b) 
demonstrate that gas drag regulates eccentricity behav- 
ior at 2-3 AU only for bodies smaller than 1 km, which 
is below our adopted threshold size d s . 

However, gas drag may have detrimental effect on cir- 
cumbinary planet formation if the turbulence (Meschiari 
2012b) or the density waves driven by the binary in the 
disk (Marzari et al. 2008) can strongly excite planetesi- 
mal eccentricities at large r/a b . Drag-induced inspiral of 
solids is another possible obstacle for planet formation. 
We leave detailed exploration of these issues to future 
study. 
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APPENDIX 

DISK-INDUCED PRECESSION RATES. 

Potential Ud of an axisymmetric disk with a power-law surface density profile in the form (|16|) . truncated at the 
inner radius r,- n (i.e. £(r) = for r < rj n ) is given by the following expression far from the inner edge of the disk, 
r > r m (Ward 1981): 

^ = | (2i + /^ p + 2) ~«. (AD 

where Ai = [(2Z) !/ (2 2Z (Z !) 2 )] 2 , and the numerical value of Kd is for p — 3/2. Substituting this expression into equation 
((U), assuming p = 3/2 and using equation (flT)) . we arrive at the result <jTUJ) . 

To calculate apsidal precession of the binary we use the expression for Ud inside the inner edge (r < fj n ) of a sharply 
truncated disk (Ward 1981): 

/ \ 00 A 21 

U d (r) = -2nGZ in r in <t> ( — ) , 0(z) = V - ' Z (A2) 

\ r in J ; „ Zl -\- p 1 

Taking into account that both stars move in this potential one can show using equation Q that binary precesses at 
the rate (HO]) with 

^) = f>f&^j + (1 - n) 21 ] z 2 ^\ (A3) 
i=i ZL+P 1 

where equation (TT71) was used. We numerically found that 4> w 0.3 — 0.8 for a^jrin = 0.3 — 0.5 and \i = — 0.5. For 
simplicity, in this study we simply set (f> w 0.5. 



